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Tensor network decompositions offer an efficient description of certain many-body states of a lattice system 
and are the basis of a wealth of numerical simulation algorithms. In a recent paper [arXiv:0907.2994vl] we 
discussed how to incorporate a global internal symmetry, given by a compact, completely reducible group Q, 
into tensor network decompositions and algorithms. Here we specialize to the case of Abelian groups and, 
for concreteness, to a [7(1 ) symmetry, often associated with particle number conservation. We consider tensor 
networks made of tensors that are invariant (or covariant) under the symmetry, and explain how to decompose 
and manipulate such tensors in order to exploit their symmetry. In numerical calculations, the use of U{\) 
symmetric tensors allows selection of a specific number of particles, ensures the exact preservation of particle 
number, and significantly reduces computational costs. We illustrate all these points in the context of the multi- 
scale entanglement renormalization ansatz. 

PACS numbers: 03.67.-a, 03.65.Ud, 03.67 .Hk 



I. INTRODUCTION 

Tensor networks are becoming increasingly popular as a 
tool to represent wave-functions of quantum many-body sys- 
tems. Their success is based on the ability to efficiently de- 
scribe the ground state of a broad class of local Hamiltonians 
on the lattice. Tensor network states are used both as a varia- 
tional ansatz to numerically approximate ground states and as 
a theoretical framework to characterize and classify quantum 
phases of matter. 

Examples of tensor network states for one dimensional 
systems include the matrix product state^l (MPS), which 
results naturally from both Wilson's numerical renormal- 
ization groups and White's density matrix renormalization 
group 5 8 (DMRG) and is also used as a basis for simula- 
tion of time evolutionpKB] the tree tensor network^] (TTN), 
which follows from coarse-graining schemes where the spins 
are blocked hierarchically; and the multi-scale entangle- 
ment renormalization ansatiLMlD(MERA), which results from 
a renormalization group procedure known as entanglement 
renormalization. For two dimensional lattices, there are 
generalizations of these three tensor network states , nam ely 
projected entangled pair states^EQ (PEPS), 2D TTr42^and 
2D MERA , ' 34 " 40 ! respectively. As variational ansatze, PEPS 
and 2D MERA are particularly interesting since they can 
be used to address large two-dimensional lattices, including 
systems of frustrated spin d 31 ! 40 ! and interacting fermionsp2E2l 
where Monte Carlo techniques fail due to the sign problem. 

A many-body Hamiltonian H may be invariant under cer- 
tain transformations, which form a group of symmetries. 50 
The symmetry group divides the Hilbert space of the the- 
ory into symmetry sectors labeled by quantum numbers or 
conserved charges. On a lattice one can distinguish between 
space symmetries, which correspond to some permutation of 
the sites of the lattice, and internal symmetries, which act on 
the vector space of each site. An example of space symme- 
try is invariance under translations by some unit cell, which 
leads to conservation of momentum. An example of internal 
symmetry is SU(2) invariance, e.g. spin isotropy in a quantum 
spin model. An internal symmetry can in turn be global, if it 
transforms the space of each of the lattice sites according to 



the same transformation (e.g. a spin independent rotation); or 
local, if each lattice site is transformed according to a differ- 
ent transformation (e.g. a spin-dependent rotation), as it is in 
the case of gauge symmetric models. A global internal SU(2) 
symmetry gives rise to conservation of total spin. By target- 
ting a specific symmetry sector during a calculation, computa- 
tional costs can often be significantly reduced while explicitly 
preserving the symmetry. It is therefore not surprising that 
symmetries play an important role in numerical approaches. 

In tensor network approaches, the exploitation of global in- 
ternal s ymmetri es has a long history, especially in the context 
of MPS pEEHnEH Both Abelian and non-Abelian symmetries 
have been thoroughly incorporated into DMRG code and have 
been exploited to obtain computational gains. Symmetries 
have also been used in more recent proposals to simulate time 
evolution with MPS, e.g. with the time evolving block dec- 
imation (TEBD) algorithm and variations thereof, often col- 
lectively referred to as time-dependent DMRG. 

When considering symmetries, it is important to notice that 
an MPS is a trivalent tensor network. That is, in an MPS 
each tensor has at most three indices. The Clebsch-Gordan 
coefficients 50 (or coupling coefficients) of a symmetry group 
are also trivalent, and this makes incorporating the symme- 
try into a MPS by considering symmetric tensors particularly 
simple. In contrast, tensor network states with a more elabo- 
rated network of tensors, such as MERA or PEPS, consist of 
tensors having a larger number of indices. In this case a more 
general formalism is required in order to exploit the symme- 
try. As explained in Ref. 66, a generic symmetric tensor can 
be decomposed into a degeneracy part, which contains all de- 
grees of freedom not determined by symmetry, and a struc- 
tural part, which is completely determined by symmetry and 
can be further decomposed as a trivalent network of Clebsch- 
Gordan coefficients. 

The use of symmetric tensors in more complex tensor net- 
works has also been discussed in Refs. 67 68 In particular, 
Ref. [67] has shown that under convenient conditions (injec- 
tivity), a PEPS that represents a symmetric state can be rep- 
resented with symmetric tensors, generalizing similar results 
for MPS obtained in Ref. |60j Notice that these studies are not 
concerned with how to decompose symmetric tensors so as 



to computationally exploit the symmetry. On the other hand, 
exploitation of U(l) symmetry for computational gain in the 
context of PEPS was reported in Ref. |68] although no imple- 
mentation details were provided. Finally, several aspects of 
local internal symmetries in tensor networks algorithms have 
been addressed in Refs. l69if72l 

The purpose of this paper is to address, in considerable de- 
tail and at a pedagogical level, several practical aspects of 
the exploitation of global internal symmetries not covered in 
Ref. |66] For concreteness we will concentrate on the U(l) 
symmetry, but extending our results to any Abelian group is 
straightfoward. A similar analysis of non-abelian groups will 
be considered in Ref. [73] 

The paper is organized in sections as follows. Section [TT] 
contains a review of the tensor network formalism and intro- 
duces the nomenclature and diagrammatical representation of 
tensors used in the rest of the paper. It also describes a set 
V of primitives for manipulating tensor networks, consisting 
of manipulations that involves a single tensor (permutation, 
fusion and splitting of the indices of a tensor) and matrix op- 
erations (multiplication and factorization). 
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Section III reviews basic notions of representation theory of 
the Abelian group U(l). The action of the group is analysed 
first on a single system, where t/(l) symmetric states and t/(l) 
invariant operators are decomposed in a compact, canonical 
manner. This canonical form allows us to identify the degrees 
of freedom which are not constrained by the symmetry. The 
action of the group is then also analysed on the tensor product 
of two Hilbert spaces and, finally, on the tensor product of a 
finite number of spaces. 

Section|lV]explains how to incorporate the t/(l) symmetry 
into a generic tensor network algorithm, by considering U(l) 
invariant tensors in a canonical form, and by adapting the set 
V of primitives for manipulating tensor networks. These in- 
clude the multiplication of two t/(l) invariant matrices in their 
canonical form, which is at the core of the computational sav- 
ings obtained by exploiting the symmetry in tensor network 
algorithms. 

Section [V] illustrates the practical exploitation of the f/(l) 
symmetry in a tensor network algorithm by presenting MERA 
calculations of the ground state and low energy states of two 
quantum spin chain models. Section|Vl]contain some conclu- 
sions. 

The canonical form offers a more compact description of 
U(l) invariant tensors, and leads to faster matrix multipli- 
cations and factorizations. However, there is also an addi- 
tional cost associated with mantaining an invariant tensor in 
its canonical form while reshaping (fusing and/or splitting) its 
indices. In some situations, this cost may offset the benefits 
of using the canonical form. In the appendix we discuss a 
scheme to lower this additional cost in tensor network algo- 
rithms that are based on iterating a sequence of transforma- 
tions. This is achieved by identifying, in the manipulation 
of a tensor, operations which only depend on the symmetry. 
Such operations can be precomputed once at the beginning of 
a simulation. Their result, stored in memory, can be re-used at 
each iteration of the simulation. The appendix describes two 
such specific precomputation schemes. 



(«) 



O rO A 



FIG. 1: (i) Graphical representation of a tensor f of rank k and com- 
ponents Tj^ 2 ... jr The tensor is represented by a shape (circle) with 
k emerging lines corresponding to the k indices i\,i2, • • ■ , h- Notice 
that the indices emerge in counterclockwise order, (ii) Graphical rep- 
resentation of tensors with rank k = 0, 1 and 2, corresponding to a 
complex number c e C, a vector |v) 6 C' 1 ' and a matrix M e C |,l|x| ' 21 , 
respectively. 



II. REVIEW: TENSOR NETWORK FORMALISM 

In this section we review background material concerning 
the formalism of tensor networks, without reference to sym- 
metry. We introduce basic definitions and concepts, as well 
as the nomenclature and graphical representation for tensors, 
tensor networks, and their manipulations, that will be used 
throughout the paper. 

A. Tensors 



A tensor f is a multidimensional array of complex num- 
bers Zf,i 2 ...4 e C. The rank of tensor T is the number k of 
indices. For instance, a rank-zero tensor (k — 0) is a complex 
number. Similarly, rank-one (k = 1) and rank-two (k = 2) 
tensors represent vectors and matrices, respectively. The size 
of an index i, denoted is the number of values that the in- 
dex takes, i e {1,2, ■• • , |/|). The size of a tensor T, denoted 
|7*|, is the number of complex numbers it contains, namely 
|f| = |/i|x|i 2 |X"-x|i k |. 

It is convenient to use a graphical representation of tensors, 
as introduced in Fig. [I] where a tensor f is depicted as a cir- 
cle (more generally some shape, e.g. a square) and each of 
its indices is represented by a line emerging from it. In order 
to specify which index corresponds to which emerging line, 
we follow the prescription that the lines corresponding to in- 
dices {z"i , iz, ■ ■ ■ , ik) emerge in counterclockwise order. Unless 
stated otherwise, the first index will correspond to the line 
emerging at nine o'clock (or the first line encoutered while 
proceeding counterclockwise from nine o'clock). 

Two elementary ways in which a tensor T can be trans- 
formed are by permuting and reshaping its indices. A permu- 
tation of indices corresponds to creating a new tensor T' from 
T by simply changing the order in which the indices appear, 
e.g. 



en 



acb — -* abc 



(i) 



On the other hand, a tensor t can be reshaped into a new 
tensor T' by 'fusing' and/or 'splitting' some of its indices. 
For instance, in 



(f')ad = T abc , d = bxc, 



(2) 



tensor t' is obtained from tensor t by fusing indices b e 
{!,■■• , \b\) and c e {!,•■■, \c\) together into a single index 
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FIG. 2: Transformations of a tensor: (z) Permutation of indices b and 
c. (ii) Fusion of indices b and c into d = b x c; splitting of index 
d = bxc into £> and c. 
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FIG. 3: (i) Graphical representation of the matrix multiplication of 
two matrices R and § into a new matrix t Q (ii) Graphical represen- 
tation of an example of the contraction of two tensors R and S into a 
new tensor t {5}. 

of of size \d\ - \b\ ■ \d\ that runs over all pair of values of b and 
c, i.e. d e {(1, 1), (1, 2), ■ • ■ , (|fc|, |c| - 1), (|i|, |c|)}, whereas in 

fabc = {f')ad, d = bXC, (3) 

tensor T is recovered from T' by splitting index c/ of T' back 
into indices b and c. The permutation and reshaping of the 
indices of a tensor have a straightforward graphical represen- 
tation; see Fig. [2] 



Again the multiplication of two tensors can be graphically rep- 
resented by connecting together the lines corresponding to in- 
dices that are being contracted (indices b and c in Eq. B); see 
Fig. git). 

The multiplication of two tensors can be broken down into 
a sequence of elementary steps to transform the tensors into 
matrices, multiply the matrices, and transform the resulting 
matrix into a tensor. Next we describe these steps for the con- 
traction given in Eq.|5] They are illustrated in Fig. [4] 

1 . Permute the indices of tensor R in such a way that the 
indices to be contracted, b and c, appear in the last posi- 
tions and in a given order, e.g. be; similarly, permute the 
indices of S so that the indices to be contracted, again b 
and c, appear in the first positions and in the same order 
be: 

(R')ad be = Rabcd 

(S')bc fh ~ S c fbh (6) 

2. Reshape tensor R' into a matrix R" by fusing into a sin- 
gle index u all the indices that are not going to be con- 
tracted, u — a x d, and into a single index y all indices 
to be contracted, y = bxc. Similarly, reshape tensor S' 
into a matrix S " with indices y - bxc and w — f xh, 

(R")uy = (R')adbc 

(S")yw = (S')bcfh- (7) 

3. Multiply matrices R" and S " to obtain a matrix T", with 
components 

(f"U = J](R")uy (8) 



B. Multiplication of two tensors 

Given two matrices R and S with components R a t, and S he, 
we can multiply them together to obtain a new matrix T,T = 
R ■ S with components 



b 



ab'J be, 



(4) 



by summing over or contracting index b. The multiplication 
of matrices R and S is represented graphically by connecting 
together the emerging lines of R and S corresponding to the 
contracted index, as shown in Fig. [3ji). 

Matrix multiplication can be generalized to tensors. For 
instance, given tensors R and S with components R„bcd and 
S c fbh> we can define a tensor T with components Tj w fd given 
by 



1 hafd 



/ \ RabcdS 
be 



cfbh- 



(5) 



4. Reshape matrix f" into a tensor f' by splitting indices 
u — axd and w - f xh, 



(T')adfh - (T")uw 



(9) 



5. Permute the indices of f' into the order in which they 
appear in T, 



Thafd - (T')adfh- 



(10) 



We note that breaking down a multiplication of two tensors 
into elementary steps is not necessary - one can simply imple- 
ment the contraction of Eq. [5] as a single process. However, 
it is often more convenient to compose the above elementary 
steps since, for instance, in this way one can use existing lin- 
ear algebra libraries for matrix multiplication. In addition, it 
can be seen that the leading computational cost in multiply- 
ing two large tensors is not changed when decomposing the 



contraction in the above steps. In Sec. IV I this subject will be 
discussed in more detail for £/(!) invariant tensors. 
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FIG. 4: Graphical representations of the five elementary steps 1-5 
into which one can decompose the contraction of the tensors of Eq.|5] 




FIG. 5: (i) Factorization of a matrix f according to a singular value 
decomposition l | 1 I [ i. (ii) Factorization of a rank-4 tensor f according 
to one of several possible singular value decompositions. 

C. Factorization of a tensor 

A matrix T can be factorized into the product of two (or 
more) matrices in one of several canonical forms. For in- 
stance, the singular value decomposition 

fab = ^ U ac S cd Vdb = 2j U ac s c V ch (11) 

c,d c 

factorizes T into the product of two unitary matrices U and V, 
and a diagonal matrix S with non-negative diagonal elements 
s c = S cc known as the singular values of T, see Fig. Hi). On 
the other hand, the eigenvalue or spectral decomposition of a 
square matrix T is of the form 

tab = Yj M ac D cd {M' l )db = Yj MacMM' l )cb (12) 

c,d e 

where M is an invertible matrix whose columns encode the 
eigenvectors \A C ) of T, 

t\A c ) = A C \A C ), (13) 

M~' is the inverse of M, and D is a diagonal matrix, with the 
eigenvalues A c = D cc on its diagonal. Other useful factoriza- 
tions include the LU decomposition, the QR decomposition, 




FIG. 6: (i) Example of a tensor network TV. (ii) Tensor f of which 
the tensor network TV could be a representation, (iii) Tensor T can 
be obtained from TV through a sequence of contractions of pairs of 
tensors. Shading indicates the two tensors to be multiplied together 
at each step. 



etc. We refer to any such decomposition generically as a ma- 
trix factorization. 

A tensor T with more than two indices can be converted 
into a matrix in several ways, by specifying how two join its 
indices into two subsets. After specifying how tensor T is to 
be regarded as a matrix, we can factorize T according to any 
of the above matrix factorizations, as illustrated in Fig. |5jii) 
for a singular value decomposition. This requires first per- 
muting and reshaping the indices of T to form a matrix, then 
decomposing the later, and finally restoring the open indices 
of the resulting matrices into their original form by undoing 
the reshapes and permutations. 



D. Tensor networks and their manipulation 

A tensor network TV is a set of tensors whose indices are 
connected according to a network pattern, e.g. Fig. [6] 

Given a tensor network TV, a single tensor T can be obtained 
by contracting all the indices that connect the tensors in TV. 
Here, the indices of tensor T correspond to the open indices 
of the tensor network TV. We then say that the network TV is 
a tensor network decomposition of T. One way to obtain T 
from TV is through a sequence of contractions involving two 
tensors at a time, Fig. [6] 

From a tensor network decomposition N for a tensor T, an- 
other tensor network decomposition for the same tensor T can 
be obtained in many ways. One possibility is to replace two 
tensors in TV with the tensor resulting from contracting them 
together, as is done in each step of Fig. |6|ii). Another way 
is to replace a tensor in N with a decomposion of that ten- 
sor (e.g. with a singular value decomposition). In this paper, 
we will be concerned with manipulations of a tensor network 
that, as in the case of multiplying two tensors or decomposing 
a tensor, can be broken down into a sequence of operations 
from the following list: 



1. Permutation of the indices of a tensor, Eq.[T] 

2. Reshape of the indices of a tensor, Eqs.|2]|3] 

3. Multiplication of two matrices, Eq.|4] 

4. Decomposition of a matrix (e.g. singular value decom- 
position (jTTJ or spectral decomposition ( 12 1). 

These operations constitute a set P of primitive operations 
for tensor network manipulations (or, at least, for the type of 
manipulations we will be concerned with). 



In Section IV we will discuss how this set P of primitive 
operations can be generalized to tensors that are symmetric 
under the action of the group t/(l). 
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E. Tensor network states for quantum many-body systems 

As mentioned in the introduction, tensor networks are used 
as a means to represent the wave-function of certain quantum 
many-body systems on a lattice. Let us consider a lattice X. 
made of L sites, each described by a complex vector space 
V of dimension d. A generic pure state |*P) e V 0L of £, can 
always be expanded as 



l*> = X 



JZ1>|Z2>-"I*L>, 



(14) 



where i s = 1, ■• ■ ,d labels a basis \i s ) of V for site s e £. 
Tensor W, with components y Vi l j 2 ...j L , contains d L complex co- 
efficients. This is a number that grows exponentially with the 
size L of the lattice. Thus, the representation of a generic pure 
state |*P) € V® L is inefficient. However, it turns out that an 
efficient representation of certain pure states can be obtained 
by expressing tensor *F in terms of a tensor network. 

Fig. [7] shows several popular tensor network decomposi- 
tions used to approximately describe the ground states of lo- 
cal Hamiltonians H of lattice models in one or two spatial 
dimensions. The open indices of each of these tensor net- 
works correspond to the indices i\, ii, ■ ■ ■ Jl of tensor T*. No- 
tice that all the tensor networks of Fig.|7]contain O(L) tensors. 
If p is the rank of the tensors in one of these tensor networks, 
and^ is the size of their indices, then the tensor network de- 
pends on 0{Lx p ) complex coefficients. For a fixed value of 
X this number grows linearly in L, and not exponentially. It 
therefore does indeed offer an efficient description of the pure 
state |*P) G V 8i that it represents. Of course only a subset 
of pure states can be decomposed in this way. Such states, 
often referred to as tensor network states, are used as vari- 
ational ansatze, with the 0(Lx p ) complex coefficients as the 
variational parameters. 

Given a tensor network state, a variety of algorithms (see 
e.g. Refs. |4]|4"9j> are used for tasks such as: (J) computation of 
the expectation value ( v F|o| v F) of a local observable o, (if) op- 
timization of the variational parameters so as to minimize the 
expectation value of the energy (*F|i7| v I'), or (z'z'z) simulation of 
time evolution, e.g. e~ lHf |*I'). These tasks are accomplished 
by manipulating tensor networks. 



FIG. 7: Examples of tensor network states for ID systems: (i) matrix 
product state (MPS), (it) tree tensor network (TTN), (Hi) multi-scale 
entanglement renormalization ansatz (MERA). Examples of tensor 
network states for 2D systems: (iv) projected entangled-pair state 
PEPS, (v) 2D TTN. (2D MERA not depicted). 



On most occasions, all required manipulations can be re- 
duced to a sequence of primitive operations in the set P intro- 



duced in Sec. II D Thus, in order to adapt the tensor network 
algorithms of e.g. Refs. 4 49 to the presence of a symmetry, 
we only need to modify the set P of primitive tensor network 



operations. This will be done in Sec. IV 



F. Tensors as linear maps 

A tensor can be used to define a linear map between vec- 
tor spaces in the following way. First, notice that an index i 
can be used to label a basis {\i}} of a complex vector space 
y[d s f dimension |z'|. On the other hand, given a tensor 
T of rank k, we can attach a direction 'in' or 'out' to each in- 
dex i\ , ii, ■ ■ ■ , I*. This direction divides the indices of T into a 
subset / of incoming indices and the subset O of outgoing in- 
dices. We can then build input and output vector spaces given 
by the tensor product of the spaces of incoming and outgoing 
indices, 



';] y[°ut] _ 



(15) 



i,€l 



i,£0 



and use tensor f to define a linear map between V [ln] and 
y[out] p or ms t ancei if a rank-3 tensor T a b c has one incom- 
ing index eel and two outgoing indices a, b e O, then it 
defines a linear map t : V [c] -> V [fll ® V w given by 



2 t abc \a)\b){c\ 



(16) 



Graphically, we denote the direction of an index by means of 
an arrow; see Fig.|8ji). 

By decorating the lines of a tensor network N with arrows 
(Fig. [8Tii)), this can be regarded as a composition of linear 
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FIG. 8: (i) Tensor f with one incoming index and two outgoing in- 
dices, denoted by incoming and outgoing arrows respectively i ] 
(ii) A tensor network N with directed links can be interpreted as a 
linear map between incoming and outgoing spaces (of the incoming 
and outgoing indices) obtained by composing the linear maps asso- 
ciated with each of the tensors in N. 



maps — namely one linear map for each tensor in N. While 
arrows might be of limited relevance in the absence of a sym- 
metry, they will play an important role when we consider sym- 
metric tensors since they specify how the group acts on each 
index of a given tensor. 



m. REVIEW: REPRESENTATION THEORY OF THE 
GROUP U(l) 

In this section we review basic background material con- 
cerning the representation theory of the group U(l). We first 
consider the action of U( 1 ) on a vector space V, which decom- 
poses into the direct sum of (possibly degenerate) irreducible 
representations. We then consider vectors of V that are sym- 
metric (invariant or covariant) under the action of t/(l), as 
well as linear operators that are t/(l) invariant. Then we con- 
sider the action of U(l) on the tensor product of two vector 
spaces, and its generalization to the tensor product of an arbi- 
trary number of vector spaces. 



A. Decomposition into direct sum of irreducible 
representations 

Let V be a finite dimensional space and let <p e [0, 2n) label 
a set of linear transformations W, 



(17) 



that are a unitary representation of the group t/(l). That is 

Wj% = %Wj = I, V tp e [0, 2tt), (18) 

%, = fy^fyy, = tfWb, v <Pu <P2 e [0, 2tt).(19) 

Then V decomposes as the direct sum of (possibly degener- 
ate) one-dimensional irreducible representations (or irreps) of 
U(l), 

(20) 

where V„ is a subspace of dimension d„, made of d n copies 
of an irrep of U(l) with charge n e Z. We say that irrep n is 



d n -fo\d degenerate and that V„ is the degeneracy space. For 
concreteness, in this paper we identify the integer charge n as 
labelling the number of particles (another frequent identifica- 
tion is with the z component of the spin, in which case semi- 
integer numbers may be considered). The representation of 
group U( \) is generated by the particle number operator h, 



h = 2^ nP n , P n = ^ \nt n )(nt n \, 



(21) 



where P„ is a projector onto the subspace V„ of particle num- 
ber n, and the vectors \nt n ), 



h\nt n ) — n\nt n ), t n = 1, 



,d n , 



(22) 



are an orthonormal basis of In terms of «, the transforma- 

% = e-' A f. (23) 



tions read 



It then follows from Eq. 22 that 

Wy«f„> = e-^\nt n ). 



(24) 



The dual basis {(nf„|) is transformed by the dual representa- 
tion of t/(l), with elements wj , as 



(25) 



Example 1: Consider a two-dimensional space V that de- 
composes as V = Vo © Vi, where the irreps n — and n — 1 
are non-degenerate (i.e. do — d\ = 1). Then the orthogonal 
vectors {\n = 0, fo = 1), \n = 1, t\ = 1)} form a basis of V. In 
column vector notation, 



= |n = 0,f = l>, 



ee \ n =l,ti = 1), 



\0j~ "" u \1 
the particle number operator h and transformation W^, read 

(I 



(26) 




1 



w = . 

f ~ 10 e~ ilp 



(27) 



Example 2: Consider a four-dimensional space V that de- 
composes as V = Vo © Vi © V2, where do = di_ — 1 and 
d\ = 2, so that now irrep n — 1 is two-fold degenerate. Let 
\\n — 1, t\ — 1), \n = 1, t\ = 2)} form a basis of Vi. In column 
vector notation, 






0) 



1 

oj 



ee \n = Q,t = l), 



hh=2), 



^0^ 
1 


oj 

(0\ 



u 



ee \n= l,ti = 1), 



(28) 



2,? 2 = 1>. (29) 



the particle number operator h and transformation read 

^ 







(0 








0^ 









1 








, w = 


h = 








1 







lo 








2; 





1 

e^ 












e-' l2 f. 



(30) 
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B. Symmetric states and operators 

In this work we are interested in states and operators that 
have a simple transformation rule under the action of U(l). A 
pure state IT*) e V is symmetric if it transforms as 



WJV) 



(31) 



The case n — corresponds to an invariant state, Wy^P) = 
which transforms trivially under U(l), whereas for n + 
the state is covariant, with \¥) being multiplied by a non- 
trivial phase e~ mip . Notice that a symmetric state I*!*} is an 
eigenstate of ft: that is, it has a well-defined particle number 
n. I*!*) can thus be expanded in terms of a basis of the relevant 
subspace V„, 

d„ 

h\¥) = n\y„), m = Yp^nWntn). (32) 



A linear operator t : 
with the generator ft , 



/„=! 



is invariant if it commutes 



[7\n]=0, (33) 
or equivalently if it commutes with the action of the group, 



w/w; = t 



G [0, 2tt). 

It follows that T decomposes as (Schur's lemma) 
f = fP>f„ 



where f„ is a d„ x d„ matrix that acts on the subspace 
Eq.m 



(34) 



(35) 



in 



Notice that the operator T in Eq. 35 transforms vectors with 
a well defined particle number « into vectors with the same 
particle number. That is, i/(l) invariant operators conserve 
particle number. 

Example 1 revisited: In Example 1 above, symmetric 
vectors must be proportional to either \n = 0, to = 1) or 
\n = \,t\ = 1). An invariant operator T = To ffi T\ is of 
the form 



a 
aj 



aro, a\ G 



(36) 



Example 2 revisited: In Example 2 above, a symmetric 
vector \W) must be of the form 



(37) 









'0) 




(0) 


PF> = 





loj 




a\ 
Pi 

Loj 


, or \Y) = 





fill, 



where ao, a\,P\, ai £ 
7*2 is of the form 



An invariant operator t = Tq © f i © 











^ 





a j 


Pi 








n 


Si 

















r = n I' ^ n (38) 



where fi corresponds to the 2x2 central block and 
ao,ai,Pi,yi,di,a2 G C. 

The above examples illustrate that the symmetry imposes 
constraints on vectors and operators. By using an eigenbasis 
{|«f„)} of the particle number operator ft, these constraints im- 
ply the presence of the zeros in Eqs. [36p8| Thus, a reduced 
number of complex coefficients is required in order to describe 
i/(l) symmetric vectors and operators. As we will discuss in 



Sec. IV performing manipulations on symmetric tensors can 
also result in a significant reduction in computational costs. 



C. Tensor product of two representations 



Let Y ( ' and V 1 ' be two spaces that carry representations 
of t/(l), as generated by particle number operators n (A) and 
n (B \ and let 



(39) 



be their decompositions as a direct sum of (possibly degen- 
erate) irreps. Let us also consider the action of U(l) on the 
tensor product Y (AB) = V (A) ® V (B) as generated by the total 
particle number operator 



. + . 



) n 



that is, implemented by unitary transformations 



(40) 



(41) 



The space V ( ' also decomposes as the direct sum of (pos- 
sibly degenerate) irreps, 



(42) 



Here the subspace V^f , with total particle number n^B, cor- 
responds to the direct sum of all products of subspaces V| 



(A) 



and 



ytB) 



such that + «b = n^B, 



" n A w v « B 



For each subspace V^f' 
basis \\n AB t„ AB )}, 



(43) 

we introduce a coupled 
(AB) \n AB t nAB ) = n AB \n AB t nAB ), (44) 



in Eq. 



42 



where each vector \n A Bt nAB ) corresponds to the tensor product 
\n A tn A \n B t nB ) = \n A t„ A ) ® \n B t„ B ) of a unique pair of vectors 
\n A t„ A ) and \n B t„ B ), with n A + n B - n A s. Let table T f " sc , with 
components 



refuse 



= (n A Bt„ AB \n A t n ; n B t„ B }, 



(45) 



encode this one-to-one correspondence. Notice that each 
component of Y |USL is either a zero or a one. Then 



\n AB t„ AB ) = ^ 



[« A ^;nBf« B >- (46) 



n A t„ A n B t„ B 
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For later reference, we notice that Y' 1 ™ can be decom- 
posed into two pieces. The first piece expresses a basis 
{\n A t„ A ; n B t„ B )} of Y (AB) in terms of the basis {\n A t„ A )} of V (A) 
and the basis {\riBt„ B )} of V (B) . This assignment occurs as in 
the absence of the symmetry, where one creates a composed 
index d — b x c by running fast over index c, as for example 
in Eq. [2] The second piece is a permutation of basis elements 
that reorganizes them according to their total particle number 
n AB . Finally, the product basis can be expressed in terms of 
the coupled basis 



induces a decomposition 



\n,\ 



with 



•„ AB ^n A t„ A ,n B t„ B 



\n AB tn AB ), (47) 



= T 



n A t„ A ,n B t„ B — >n/lfl'n /1B ' 



(48) 



t^t split 

" A Bt,i AB ->n A t„ A ,n B t„ B 



Example 3: Consider the case where both V ( ' and 
correspond to the space of Example 1, that is V (A) = V^'ffiV^ 
i \W Y[ B \ where V[, A) , Vf\ Y (B \ and vf > all 



and 

have dimension one. Then 
Example 2, namely 



corresponds to the space in 



where 



AAB) 




? (AB) ( 



r(B> 




(49) 



(50) 



V^sY^V™. (52) 
The coupled basis {InAB^s)) reads, 



|n A B=0,? = l> = |ba = 0,/b = 1>«|b b = 0,* = 1> (53) 

\n AB = l,*i = 1> = = 0,t = l)®\n B = l,h = 1> (54) 

\n AB = l,fi = 2) = \n A = l,h = \)®\n B =0,t = 1) (55) 

\n AB =2,t 2 = 1> = I«a = Mi = l)®|n B = l,h = 1), (56) 

where we emphasize that the degeneracy index t„ w takes two 
possible values for n AB = 1, i.e. t\ e {1,2}, since there are 
two states \n A t„ A )® \n B t„ B ) with n A + n B — 1. The components 
^n A t A ,n B t B ^n AB t AB °f tne t ensor T' USL ' that encodes this change of 
basis all zero except for 

1 01,01^01 _ 1 01, 11^11 _ 1 11,01— »12 _ 1 11,1 1— »21 _ i - 



D. Lattice models with 1/(1) symmetry 

The action of U(l) on the three-fold tensor product 

Y (ABC) = V (A) ®V (B) ®V (C) , (57) 
as generated by the total particle number operator 



(ABC) _ "(A) 



= h (A) ® I ® I + I ® h (B) ® I + I® 1® n (C) , (58) 



(59) 



in terms of irreps V^JJ? which we can now relate to V„ A , 



-aB) 



and 



AC) 



„, c.i^. ,, . For example, we can first consider the product 

<? = V<>V^ and then the product V^ C) = V<£f ®v£>, 
and use two tables T'" sc to relate at each step the coupled basis 
with the product basis, as discussed in the previous section. 
Similarly we could consider the action of t/(l) on four tensor 
products, and so on. 

In particular we will be interested in a lattice _£ made of L 
sites with vector space V 8i , where for simplicity we assumed 
that each site s e X. is described by the same finite dimen- 
sional vector space V (see Sec. [HE}. Given a particle number 
operator h defined on each site, we can consider the action of 
U( 1 ) generated by the total particle number operator 



S=] 



which corresponds to unitary transformations 
The tensor product space V 8i decomposes as 



(60) 



(61) 



(62) 



and we denote by fliVfjv)} the particle number basis in V® L . 

We say that a lattice model is t/(l) symmetric if its Hamil- 
tonian H : V — > V commutes with the action of the group. 
That is, 

[#,#] = (63) 

or equivalently, 

(%f L 6(Wjf L = 6 \/<pe[0,2n). (64) 

One example of a t/(l) symmetric model is the Hardcore 
Bose Hubbard Model, with Hamiltonian 

L L 

Hhcbh = 2^ ( a J»s+i + "AI+i + 7» A+i) - ^ n s , (65) 

1=1 s=\ 

where we consider periodic boundary conditions (by identify- 
ing sites L + 1 and 1) and a\,a s are hardcore bosonic creation 
and annihilation operators respectively. In terms of the basis 
introduced in Example 1 , these operators are defined as 



1 





1 



To see that Hhcbh commutes with the action of the group we 
first observe that for two sites 



|ajfl2 + a} 2 a\ , hi + #2] = 0, 



(66) 
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from which it readily follows that [i/#cB»>-W] = 0. 

Notice that the chemical potential term ~[i £ s h s = -fiN 
also commutes with the rest of the Hamiltonian. The ground 
state ITS) of Hhcbh m a particular subspace or particle 
number sector can be turned into the absolute ground state by 
tuning the chemical potential fi. This fact can be used to find 
the ground state \ x ¥'^) of any particle number sector through 
an algorithm that can only minimize the expectation value of 
Hhcbh- However, we will later see that the use of symmetric 
tensors in the context of tensor network states will allow us to 
directly minimize the expectation value of Hhcbh in a given 
particle number sector by restricting the search to states 



t N =\ 

with the desired particle number N. 
Finally, by making the identifications 



0% + ICTy 



2 2 
where <x v , <x v , & z are the Pauli matrices 



1 

1 0/' 



(Ty = 



-i 

1 



1 
-1 



(68) 



one can map Hhcbh to the spin-^ XXZ quantum spin chain 

H XXZ . J] (<5f ^ +1) + dfd?» + MV>d?») , (69) 



where we have ignored terms proportional to N and A = y/4. 
In particular, for A = we obtain the quantum XX spin chain 



(70) 



1=1 

and for y — 1 , the quantum Heisenberg spin chain 

Hxxx = J] + dfd?* + <^ s+1) ) • (71) 



In Sec.|V] the quantum spin models (70 1 and (71 1 will be used 
to benchmark the performance increase resulting from use of 
symmetries in tensor networks algorithms. 



IV. TENSOR NETWORKS WITH U(l) SYMMETRY 

In this section we consider £7(1) symmetric tensors and ten- 
sor networks. We explain how to decompose C7(l) symmetric 
tensors in a compact, canonical form that exploits their sym- 
metry. We then discuss how to adapt the set f of primitives for 
tensor network manipulations in order to work in this form. 
We also analyse how working in the canonical form affects 
computational costs. 



A. U(l) symmetric tensors 



L et f be a rank-fc tensor with components f,,,,...,,. As in 



Sec. 



II F 



tor spaces 



we regard tensor T as a linear map between the vec- 

out] 



and 



15 i. This implies that each index 



is either an incoming or outgoing index. On each space V M , 
associated with index //, we introduce a particle number oper- 
ator n (/) that generates a unitary representation of {7(1) given 
by matrices = e~ w " >lf , tp e [0, 2n). In the following, we 
use wf * to denote the complex conjugate of wf. 



Let us consider the action of {7(1) on the space 



(67) given by 



V 1 " 1 8 V ,i:l 



where 



yd) _ 



w, 



w, 



(0 



if U e I, 
if k e O, 



(72) 



(73) 



(74) 



That is, Xr acts differently depending on whether index i; is 
an incoming or outgoing index of T. We then say that tensor 
T, with components r (1 ,v.., t , is {7(1) invariant if it is invariant 
under the transformation of Eq. [73] 

Z (^) jlh (^) 4fc "-( 4 ? , ) 5t ^ = ^M. ™ 

h,h,- ,k 

for all tp e [0, 2n). This is depicted in Fig. [9] 

Example 4: A {7(1) invariant vector \W) — that is, a vec- 
tor with h^) = and components (*P„=o)f in the subspace 
V„=o corresponding to vanishing particle number n — (cf. 
Eq.[32> — fulfills 



(%=oV = £ (%\ B% (^=o) ( „ V tp e [0, 2n), (76) 



in accordance with Eq. 31 as shown in Fig. [9] 

Example 5: A U(l) invariant matrix f , Eq.B5] fulfills 



a,b 

= V^[0,2.), (78) 



in accordance with Eq. 34 see Fi g.M 

Example 6: Tensor T in Eq. [16] with components f a b c 
where a and b are outgoing indices and c is an incoming index, 
is U(l) invariant if 

a,b,c 
a,b,c 

for all tp e [0,2tt), see Fig. [9] 
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V¥n=o) (^=0) 

o 



\W1 



6 



(0 



" w. 



(ii) 




B. Canonical form for U(l) invariant tensors 

Let us now write a tensor T in a particle number basis on 
each factor space in Eq. 72 That is, each index i\, 12, ■ ■ ■ , ik is 
decomposed into a particle number index n and a degeneracy 
index f„, i\ = (rei, t m ), i 2 = (n 2 , f« : ), • • • , h = f« t X and 



FIG. 9: (i) Constraint fulfilled by a 1/(1) invariant vector. The only 
allowed particle number on the single index is n = 0. (ii) Constraint 
fulfilled by a 1/(1) invariant matrix. It follows from Schur's lemma 
that the matrix is block-diagonal in particle number, (iii) Constraint 
fulfilled by a rank-three tensor with one incoming index and two out- 
going indices. 



(83) 



hp 



|<F) 

o 



-mcp 



(0 



= Oof J 

(b) 



Here, for each set of particle numbers «i,«2, ■ ■ ■ ,nu we re- 
gard f nxm -n k as a tensor with components (f n]n2 ...„ t ) 

Let Ni n and N out denote the sum of particle numbers corre- 
sponding to incoming and outgoing indices, 



Nh 



= ^ re/, iV 0Ut = ^ re/. 



(84) 



n<eO 



FIG. 10: (i) f/(l) covariant vector Q, with some non-vanishing parti- 
cle number n t 0. Under the action of U( 1 ) on its index, the covariant 
vector Q acquires a phase l |8 1 [ 1. (ii) The U(l) covariant vector 
Q, with components Q ii , can be represented by a U(l) invariant ma- 
trix T with components = Q ti , where (' is a trivial index (|i| = 1) 
with charge 11. 



Further, we say that a tensor Q, with components 2, ll2 ..., t , is 



The condition for a non-vanishing tensor of the form f„ inn ... nt 
to be invariant under t/(l), Eq. [73] is simply that the sum of 
incoming particle numbers equals the sum of outgoing particle 
numbers. Therefore, a t/(l) invariant tensor T satisfies 



(85) 



(We use the direct sum symbol (J) to denote that the different 
tensors T, hm ...„ t are supported on orthonormal subspaces of 



the tensor product space of Eq. 72 ) In components, the above 
expression reads, 



t/(l) covariant if under the transformation of Eq. 73 it simply 
aquires a non-trivial phase e~ m,p , 



Ti\h-i k — {Tn\n 2 -nk) t ( <^A 



(86) 



/ \ ffip^fi {^^\i> '"{^v^fi 0iik-ik ~ e m>P Q'Vi-'k' Here, 5N in ,N m implements particle number conservation: if 



for all <p e [0, 2tt). 

Example 7: A U(l) covariant vector \W) — that is, one 
which satisfies n| v F) = relT*) for some n + 0, and has 
nonzero components (¥„);„ on ly m me relevant subspace V„ 
(cf. Eq.[32t— fulfills 



Njn + N out , then all components of T nim ... nk must vanish. This 
generalizes the block structure of U(l) invariant matrices in 



Eq. 35 (where t nn is denoted T„) to tensors of arbitrary rank 



k. The canonical decomposition in Eq. 85 is important, in that 
it allows us to identify the degrees of freedom of tensor T that 
are not determined by the symmetry. Expressing tensor T in 
terms of the tensors T ni „ T .. nk with N m = N out ensures that we 



Yj(%) , QVnh = e'^QVnX' Vy>e[0,2^), (81) store T in the most compact possible way. 



in accordance with Eq. |3T| 

Notice that we can describe the rank-£ covariant tensor Q 
above by a rank-(fe +1) invariant tensor T with components 



k i - Qi\ 



(82) 



Notice that the canonical form of Eq. 85 is a particular case 
of the canonical form presented in Eq. 15 of Ref. 66 for more 
general (possibly non-Abelian) symmetry groups. There, a 
symmetric tensor was decomposed into degeneracy tensors 
(analogous to tensors t ni „ 2 ...„ t in Eq. [85} an d structural tensors 
(generalizing the term 6N in ,N oul in Eq. 85 1 which can in gen- 



This is built from Q by just adding an extra incoming index i, 
where index i has fixed particle number n and no degeneracy 
(i.e., i is associated to a trivial space V w = C). We refer to 
both invariant and covariant tensors as symmetric tensors. By 
using the above construction, in this work we will represent all 
t/(l) symmetric tensors by means of t/(l) invariant tensors. In 
particular, we represent the non-trivial components CP„)«„ of 
the covariant vector J¥ n ) in Eqs. 3Tp2 



eral be expanded as a trivalent network of Clebsch-Gordan (or 
coupling) coefficients of the symmetry group. In the case of 
non-Abelian groups, where some irreps have dimension larger 
than one, the structural tensors are highly non-trivial. How- 
ever, for the group (7(1) discussed in this paper (as for any 
other Abelian group) all irreps are one-dimensional and the 
structural tensors are always reduced to a simple expression 
such as 6ff ir „N au 



in Eq. 85 



as an invariant matrix 



T of size |f„|x 1 with components T tiil = C¥„) tn - Consequently, 
from now on, we will mostly consider only invariant tensors. 



(Nevertheless, in the appendix we 
will resort to a more elaborate decomposition of the structural 
tensors in order to further exploit the symmetry during tensor 
network manipulations of iterative algorithms.) 
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FIG. 11: A tensor network N made of (7(1) invariant tensors repre- 
sents a (7(1) invariant tensor T. This is seen by means of two equal- 
ities. The first equality is obtained by inserting resolutions of the 
identity I = on each index connecting two tensors in N. The 

second equality follows from the fact that each tensor in TV is (7(1) 
invariant. 



C. U(l) symmetric tensor networks 



In Sec. |II F| we saw that a tensor network TV where each 
line has a direction (represented with an arrow) can be inter- 
preted as a collection of linear maps composed into a single 
linear map T of which AT is a tensor network decomposition. 
By introducing a particle number operator on the vector space 
associated to each line of TV, we can define a unitary represen- 
tation of £7(1) on each index of each tensor in TV Then we say 
that TV is a U(l) invariant tensor network if all its tensors are 
t/(l) invariant. Notice that, by construction, if AT is a t/(l) in- 
variant tensor network, then the resulting linear map T is also 



t/(l) invariant. This is illustrated in Fig. 11 



More generally, we can consider a (7(1) symmetric tensor 
network, made of tensors that are t/(l) symmetric (that is, ei- 
ther invariant or covariant). Recall, however, that any covari- 
ant tensor can be represented as an invariant tensor by adding 
an extra index (82 1. Therefore without loss of generality we 



can restrict our attention to invariant tensor networks. 



D. Tensor network states and algorithms with (7(1) symmetry 



in the canonical form of Eqs. 85]|86 We next explain how to 
adapt the set P of four primitive operations for tensor network 



manipulation discussed in Sect II D namely permutation and 



reshaping of indices, matrix multiplication, and factorization. 



E. Permutation of indices 

Given a U(l) invariant tensor T expressed in the canonical 
form of Eqs. 85p6 permuting two of its indices is straight- 
foward. It is achieved by swapping the position of the two par- 
ticle numbers of T„ l „ 1 ...„ l , involved, and also the corresponding 
degeneracy indices. For instance, if the rank-3 tensor T of 



Eq. 16 is U{\) invariant and has components 



T a bc — \ Tn A n B nc) t , , "n 

N ' hi\ 'iiu hie- 



(88) 



when expressed in the particles number basis a = (ha, t„ A ), 
b = (riB,t„ B ), c = (rtc,t nc ), then tensor t' of Eq. [I] obtained 
from T by permuting the last two indices, has components 



(T )acb - \T n n n \ 5„ 

x ' 'n\ *Hr- 'no 



(89) 



Notice that since we only need to permute the components 
of those T„ A n B n c sucn that tia + n% = nc, implementing the 
permutation of indices requires les computational time than a 
regular index permutation. This is shown in Fig. 12 corre- 



sponding to a permutation of indices using MATLAB. 



F. Reshaping of indices 

The indices of a £7(1) invariant tensor can be reshaped 
(fused or split) in a similar manner to those of a regular ten- 
sor. However, maintaining the convenient canonical form of 
Eqs. [85|86| requires additional steps. Two adjacent indices 
can be fused together using the table T fusc of Eq. 45 which is 
a sparse tensor made of ones and zeros. Similarly an index 
can be split into two adjacent indices by using its inverse, the 



As discussed in Sec. ITTeI a tensor network TV can be used sparse tensor Y* plu of Eq. 48 



to describe certain pure states [W) e V of a lattice £.. If TV is 
a f7(l) symmetric tensor network then it will describe a pure 
state I*! 7 } that has a well-defined total particle number TV. That 
is, a t7(l) symmetric pure state 



Example 8 : Let us consider again the rank-3 tensor T of 
Eq. 16 with components given by Eq. 88 where a and b are 
outgoing indices and c is an incoming index. We can fuse 
outgoing index b and incoming index c into an (e.g. incoming) 
index d, obtaining a new tensor T' with components 



-Wtp 



,-iNip 



(87) 



In this way we can obtain a more refined version of popu- 
lar tensor network states such as MPS, TTN, MERA, PEPS, 
etc. As a variational ansatz, a symmetric tensor network state 
is more constrained than a regular tensor network state, and 
consequently it can represent less states \V) e V® L . How- 
ever, it also depends on less parameters. This implies a more 
economical description, as well as the possibility of reducing 
computational costs during its manipulation. 

The rest of this section is devoted to explaining how one 
can achieve a reduction in computational costs. This is based 
on storing and manipulating (7(1) invariant tensors expressed 



(f')ad - (f'„ AnD ) t t 5„ 



(90) 



where no = -rig + «c- [The sign in front of n# comes from 
the fact that d is an incoming index and b an outgoing index.] 
The components of T' are in one-to-one correspondence with 
those of T and follow from the transformation 

yn A n D ). . - /, \jn A nBnc). . , ^n B t nB ,n c t, 

'" A '"D "» A '"b'"C " 



> c -"'d;, d 



(91) 

where only the case n A -n D needs to be considered. To com- 
plete the example, let us assume that index a is described by 
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1. Then 



and 



the vector space V ( ' = Vo©Vi©V2 with degeneracies do - 1, 
d\ — 2 and dz — 1; index b is described by a vector space 
V (B) = V_i V without degeneracies, that is d-\ = do = 1; 
and index c is described by a vector space V (C) = Vo © Vi also 
without degeneracies, d-\ = do 
Eq.[9T] amounts to 

( f oo)„ - 

( f ;0 12 - 

( f n) 21 - 
( f n) 22 - 



f()oo) in , 

fl01 )m' 
*ioi) 211 . 

fl - 10 )nr 
fi-io) 211 ■ 



where we notice that tensor T' is a matrix as in Eq. 38 Sim- 
ilarly, we can split incoming index d of tensor T' back into 
outgoing index b and incoming index c of tensor T according 
to 



\Tn A n B n c ) , , ~ /_j \ T n A n D ) ^n D t„ D - 

s "»A r »B r "C ~" ' '"A T "D ° 



>n B t„ B ,n c t„ c 



(92) 



which, again, is non- trivial only for -rig + nc -no and ua + 
n B = n c . 

This example illustrates that fusing and splitting indices 
while maintaining the canonical form of Eqs. 85p6 requires 
more work than reshaping regular indices. Indeed, after tak- 
ing indices b and c into d — b x c by listing all pairs of values 
bx c, we still need to reorganize the resulting basis elements 
according to their particle number no- Although this can be 
done by following the simple table given by T ruse , it may add 
significantly to the overall computational cost associated with 



reshaping a tensor. For instance, Fig. 12 shows that, when us- 
ing MATLAB, fusing indices of invariant tensors can be more 
expensive than fusing indices of regular tensors. 



G. Multiplication of two matrices 

By permuting and reshaping the indices of a t/(l) invariant 
tensor, we can convert it into a U(\) invariant matrix T - 
©„„- frm'Snrf, or simply 



* without symmetry 
n with symmetry 



n 10 
E 



E «T 



r=0(d 4 ) 



I 10" 



•* — without symmetry 
□ with symmetry 



f=0(d 4 ) 



T = 



(93) 



FIG. 12: Computation times (in seconds) required to permute and 
fuse two indices of a rank-four tensor T, as a function of the size of 
the indices. All four indices of T have the same size 5d, and therefore 
the tensor contains \T\ = 5 4 d 4 coefficients. The figures compare the 
time required to perform these operations using a regular tensor and 
a (7(1) invariant tensor, where in the second case each index contains 
5 different values of the particle number n (each with degeneracy 
d) and the canonical form of Eqs. |85|86| is used. The upper figure 
shows the time required to permute two indices: For large d, exploit- 
ing the symmetry of a f/(l) invariant tensor by using the canonical 
form results in shorter computation times. The lower figure shows 
the time required to fuse two adjacent indices. In this case, main- 
taining the canonical form requires more computation time. Notice 
that in both figures the asymptotic cost scales as 0(d 4 ), or the size 
of T, since this is the number of coefficients which need to be rear- 
ranged. We note that the fixed-cost overheads associated with sym- 
metric manipulations could potentially vary substantially with choice 
of programming language, compiler, and machine architecture. The 
results given here show the performance of the authors' MATLAB 
implementation of 1/(1) symmetry. 



where t n = f m . In components, matrix t reads 

(f)ab = (tn) f t , 



(94) 



where a = (n, ?„) and b = (n, t' n ). In particular, similar to the 
discussion in Sec . II B for regular tensors, the multiplication of 
two tensors invariant under the action of U(l) can be reduced 
to the multiplication of two t/(l) invariant matrices. 



Let R and S be two t/(l) invariant matrices, with canonical 
forms 



tf = 0tf„, 5=05„. (95) 

n n 

Their product f = R S , Eq.[4| is then another matrix t which 
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is also block diagonal, 



T = 



(96) 



such that each block f „ is obtained by multiplying the corre- 
sponding blocks R„ and S„, 



Tn — R n ' S n . 



(97) 



Eqs. 93 and 97 make evident the potential reduction of com- 
putational costs that can be achieved by manipulating t/(l) in- 
variant matrices in their canonical form. First, a reduction in 
memory space follows from only having to store the diagonal 
blocks in Eq. 93 Second, a reduction in computational time 



is implied by just having to multiply blocks in Eq. [97] This is 
illustrated in the following example 

Example 9 : Consider a t/(l) invariant matrix T which is 
a linear map in a space V that decomposes into q irreps V„, 
each of which has the same degeneracy d„ = d. That is, T is 
a square matrix of dimensions dq x dq, and with the block- 
diagonal form of Eq. 93 Since there are q blocks T„ and each 
block has size d X d, the 1/(1) invariant matrix T contains qd 2 
coefficients. For comparison, a regular matrix of the same size 
contains q 2 d 2 coefficients, a number greater by a factor of q. 

Let us now consider multiplying two such matrices. We use 
an algorithm that requires 0{P) computational time to multi- 
ply two matrices of size I x I. The cost of performing q multi- 
scales as 0(qd y ). In con- 
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plications of d x d blocks in Eq. 
trast the cost of mutiplying two regular matrices of the same 
size scales as 0(q 3 d 3 ), requiring q 2 times more computational 
time. 
Fig. 



13 shows a comparison of computation times when 



multiplying two matrices with MATLAB, for both 1/(1) sym- 
metric and regular matrices. 



H. Factorization of a matrix 



The factorization of a 1/(1) invariant matrix T, Eq. 93 
can also benefit from the block-diagonal structure. Consider, 
for instance, the singular value decomposition T = US V of 
Eq 



1 1 In this case we can obtain the matrices 



u = Vr)u„ s = ms„ v = tfiv„ (98) 



ra 10" 



* ' • without symmetry 
o with symmetry 



□ □ - § 



t=0{<f) 



> 



E 10- 



* without symmetry 
o with symmetry 



f=0(cf) 



FIG. 13: Computation times (in seconds) required to multiply two 
matrices (upper panel) and to perform a singular value decomposi- 
tion (lower panel), as a function of the size of the indices. Matrices 
of size 5d x 5d are considered. The figures compare the time re- 
quired to perform these operations using regular matrices and 1/(1) 
invariant matrices, where for the 1/(1) matrices each index contains 
5 different values of the particle number n, each with degeneracy 
d, and the canonical form of Eqs. |93|94] is used. That is, each ma- 
trix decomposes into 5 blocks of size d x d. For large d, exploiting 
the block diagonal form of 1/(1) invariant matrices results in shorter 
computation time for both multiplication and singular value decom- 
position. The asymptotic cost scales with d as 0(d 3 ), while the size 
of the matrices grows as 0(d 2 ). We note that the fixed-cost overheads 
associated with symmetric manipulations could potentially vary sub- 
stantially with choice of programming language, compiler, and ma- 
chine architecture. The results given here show the performance of 
the authors' MATLAB implementation of U(l) symmetry. 



by performing the singular value decomposition of each block 
T„ independently, 



T„ — U n S n V n . 



(99) 



The computational savings are analogous to those described 
in Example 9 above for the multiplication of matrices. Fig. 13 
also shows a comparison of computational times required to 
perform a singular value decomposition on U(l) invariant and 
regular matrices using MATLAB. 



I. Discussion 

In this section we have seen that U(l) invariant tensors can 
be written in the canonical form of Eqs. [85p6] and that this 
canonical form is of interest because it offers a compact de- 
scription in terms of only those coefficients which are not 
constrained by the symmetry. We have also seen that main- 
taining the canonical form during tensor manipulations adds 
some computational overhead when reshaping (fusing or split- 
ting) indices, but reduces computation time when permuting 
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FIG. 14: MERA for a system ofL = 2x3 2 = 18 sites, made of two 
layers of disentanglers u and isometries w and a top tensor t . 



indices (for sufficiently large tensors) and when multiplying 
or factorizing matrices (for sufficiently large matrix sizes). 

The cost of reshaping and permuting indices is proportional 
to the size \T\ of the tensors, whereas the cost of multiplying 
and factorizing matrices is a larger power of the matrix size, 
for example |r| 3 ^ 2 . The use of the canonical form when ma- 
nipulating large tensors therefore results in an overall reduc- 
tion in computation time, making it a very attractive option in 
the context of tensor network algorithms. This is exemplified 
in the next section, where we apply the MERA to study the 
ground state of quantum spin models with a f/(l) symmetry. 

On the other hand, the cost of maintaining invariant tensors 
in the canonical form becomes more relevant when dealing 
with smaller tensors. In the next section we will also see that 
in some situations, this additional cost may significantly re- 
duce, or even offset, the benefits of using the canonical form. 
In this event, and in the specific context of algorithms where 
the same tensor manipulations are iterated many times, it is 
possible to significantly decrease the additional cost by pre- 
computing the parts of the tensor manipulations that are re- 
peated on each iteration. Precomputation schemes are de- 
scribed in more detail in the Appendices. Their performance 
is illustrated in the next section. 



V. TENSOR NETWORK ALGORITHMS WITH U(l) 
SYMMETRY: A PRACTICAL EXAMPLE 

In previous sections we have described a strategy to incor- 
porate a U(l) symmetry into tensors, tensor networks, and 
their manipulations. To further illustrate how the strategy 
works in practice, in this section we consider its implemen- 
tation in the context of the multi-scale entanglement renor- 
malization ansatz, or MERA. 



A. Multi-scale entanglement renormalization ansatz 



Fig. 
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shows a MERA that represent states IT) e V® 
of a lattice X. made of L — 18 sites (see Sec. |II E[ ). Recall 
that the MERA is made of layers of isometric tensors, known 
as disentanglers u and isometries w, that implement a coarse- 
graining transformation. In this particular scheme, isometries 



map three sites into one and the coarse-graining transforma- 
tion reduces the L = 18 sites of X. into two sites using two lay- 
ers of tensors. A collection of states on these two sites is then 
encoded in a top tensor f, whose upper index a = 1, 2, • • • ,^-, op 
is used to label x** states e V 8L . 

In this section we will consider a MERA analogous to that 



of Fig. 14 but with Q layers of disentanglers and isometries, 



which we will use to describe states on a lattice X. made of 
2 X 3 G sites. We will use this variational ansatz to obtain 
an approximation to the ground state and first excited states 
of two quantum spin chains that have a global internal U(l) 
symmetry, namely the spin- 1/2 quantum XX chain of Eq. 70 
and the spin-1 /2 antiferromagnetic quantum Heisenberg chain 
of Eq. |7T| Each spin- 1/2 degree of freedom of the chain 
is described by a vector space spanned by two orthonormal 
states {| X), | t)}. Here we will represent them by the states 
{|0), |1)} corresponding to zero and one particles, as in Exam- 
ple 1 of Sec. |III A| For computational convenience, we will 
consider a lattice £ where each site contains two spins, or 
states, (J |J,>, | |T>, I 11), I TTM- Therefore each site of £ is de- 
scribed by a space V = Vo © Vi © V2, where do = di = 1 and 



d\ = 2, as in Example 2 of Sec. Ill A Thus, a lattice _£ made 
of L sites corresponds to a chain of 2L spins. In such a system, 
the total particle number N ranges from to 2L. [Equivalently, 
the z-component of the total spin S - ranges from — L to L, with 
S Z =N-L]. 



B. MERA with U(l) symmetry 

A t/(l) invariant version of the MERA, or U(l) MERA for 
short, is obtained by simply considering U(l) invariant ver- 
sions of each isometric tensors, namely the disentanglers u, 
isometries w, and top tensor i. This requires assigning a par- 
ticle number operator to each index of the MERA. Each open 
index of the first layer of disentanglers corresponds to one site 
of £.. The particle number operator on any such index is there- 
fore given by the quantum spin model under consideration. 
We can characterize the particle number operator by two vec- 
tors ft and d — a list of the different values the particle number 
takes and the degeneracy associated with each such particle 
number, respectively. In the case of the vector space V for 
each site of £, described above, ft — [0, 1, 2] and d — [1,2, 1]. 
For the open index of the tensor f at the very top the MERA, 
the assignment of charges is also straighforward. For instance, 
to find an approximation to the ground state and first seven ex- 
cited states of the quantum spin model with particle number 
N, we choose ft = [N] and d — [8]. [In particular, a vanishing 
S z corresponds to N = L.] 

For each of the remaining indices of the MERA, the as- 
signment of the pair (ft, d) needs careful consideration and 
a final choice may only be possible after numerically test- 
ing several options and selecting the one which produces the 
lowest expectation value of the energy. Table [I] shows the as- 
signment of particle numbers and degeneracies made to rep- 
resent the ground state and several excited states in a system 
of L — 2 x 3 3 = 54 sites (that is, 108 spins) with total particle 
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number N - L - 54 [or S z = 0]. Notice that at level q of 
the MERA (q — 1, 2, 3) each index effectively corresponds to 
a block of n q = 3 q sites of £. Therefore having exactly n q 
particles in a block of n q sites corresponds to a density of 1 
particle per site of X- The assigned particle numbers of Table 
[I] namely [n q - 2,n q - \,n q ,n q + l,n q + 2] for level q, then 
correspond to allowing for fluctuations of up to two particle 
with respect to the average density. The sum of corresponding 
degeneracies <j = [d n ^ 2 , 4,-1 . d„ q , d„ q+l , d nq + z ] gives the bond 
dimension^, which in the example i&x = 13- 



Level q 


Particle numbers ft 


Degeneracy d 


top 


N = 54 




3 


[25,26, 27,28,29] 


[1,3,5,3,1] 


2 


[7,8,9,10,11] 


[1,3,5,3,1] 


1 


[1,2,3,4,5] 


[1,3,5,3,1] 





[0,1,2] 


[1,2,1] 



TABLE I: Example of particle number assignment in a U(l) MERA 
for L = 54 sites (or 108 spins). The total bond dimension is x — 
1 + 3 + 5 + 3+1 = 13. 



In order to find an approximation to the ground state of ei 



ther Hxx or Hxxx in Eqs. 70|7T we set^, op = 1 and optimize 



the tensors in the MERA so as to minimize the expectation 
value 



mam 



(100) 



where \W} e V 0i is the pure state represented by the MERA 
and H is the relevant Hamiltonian. In order to find an approx- 
imation to the Xwp > 1 eigenstates of H with lowest energies, 
we optimize the tensors in the MERA so as to minimize the 
expectation value 



■* •■ Heisenberg 
o XX 



FIG. 15: Error in ground state energy AE as a function of^ for the 
XX and Heisenberg models with 2L = 108 spins and periodic bound- 
ary conditions, in the particle number sector N = L (or S z = 0). The 
error is seen to decay exponentially with^. 



X 


Degeneracy d 


no. of coefficients 


no. of coefficients 


ratio 






(regular) 


(symmetric) 






4 


[0,1,2,1,0] 


1552 


426 


3.6 


1 


8 


[0, 2,4, 2, 0] 


17216 


4714 


3.7 


1 


13 


[1,3,5,3,1] 


115501 


21969 


5.3 


1 


17 


[1,4,7,4, 1] 


335717 


68469 


5.0 


1 


21 


[1,5,9,5,1] 


779965 


166901 


4.7 


1 


30 


[2,7, 12,7,2] 


3243076 


639794 


5.1 


1 



TABLE II: Number of coefficients required to specify the tensors of 
a MERA for L = 54 as a function of the bond dimension x, which 
decomposes into a degeneracy vector 1. A comparison is made be- 
tween regular tensors and 1/(1) invariant tensors. 
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The optimization is carried out using the MERA algorithm de- 
scribed in Ref.[T8] which requires contracting tensor networks 
(by sequentially multiplying pairs of tensors) and performing 
singular value decompositions. In the present example, all of 
these operations will be performed exploiting the f/(l) sym- 
metry. 

Fig.[l5]shows the error in the ground state energy as a func- 
tion of the bond dimension^-, for assignments of degeneracies 
similar to those in Table [II] The error is seen to decay expo- 
nentially with increasing x, indicating increasingly accurate 
approximations to the ground state. 



C. Exploiting the symmetry 

We now discuss some of the advantages of using the U(l) 
MERA. 



1. Selection of particle number sector 

An important advantage of the t/(l) MERA is that it ex- 
actly preserves the t/(l) symmetry. In other words, the states 
resulting from a numerical optimization are exact eigenvec- 
tors of the total particle number operator (60 1. In addition, 
the total particle number N can be pre-selected at the onset 
of optimization by specifying it in the open index of the top 
tensor f . 

Fig. 16 shows the energy gap between the ground state of 
an XX chain with 2L spins (or L sites), for N — L particles 
(S : = 0) and two excited states. One is the first excited state 
with also N — L particles. The other is the ground state in the 
sector with N — L + 1 particles. The two energy gaps are seen 
to decay with the system size as Zr 1 . The ability to pre-select 
a given particle number N means that only two optimizations 
were required: one MERA optimization for N — L with^f top = 
2 in order to obtain an approximation to the ground state and 
first excited state of Hxx in that particle number sector; and 
one MERA optimization for N — L + 1 with ^, op = 1 in order 
to obtain an approximation to the ground state of Hxx in the 
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exact value of A L 
exact value of A 



10' 

L 



in the particle sector N - 54 (pi S : — 0) and neighboring 
particle sectors. Recall that Hxxx is actually invariant un- 
der a global internal SU{2) symmetry, of which particle num- 
ber is a U(l) subgroup. Correspondingly the spectrum is or- 
ganized according to irreps of SU(2), namely singlets (total 
spin 0), triplets (total spin 1), quintuplets (total spin 2), etc. 
Again, using the U(l) MERA, the five particle number sectors 
N = 52,53,54,55 and 56 can be addressed with independent 
computations. This implies, for instance, that in order to find 
the gap between the first and fourth singlets, we can simply 
set N - 54 and x^ — 9 on the open index of the top tensor t. 
In order to capture the fourth singlet using the regular MERA, 
we would need to consider at least x^ v = 19 (at a larger com- 
putational cost and possibly lower accuracy), since this state 
has only the 19 th lowest energy overall. 



FIG. 16: Decay of energy gaps A with system size L in the XX 
model. The upper line corresponds to the energy gap A L between 
the ground state and the first excited state in the N = L particle num- 
ber (or S z = 0) sector. The lower line corresponds to the energy gap 
Ai + i between the ground states of the N = L and N = L + 1 particle 
number sectors. 



N 



FIG. 17: Low energy spectrum of Hxxx with L = 54 sites (=108 
spins). Depicted states have spins of zero (x), one (+), or two (o), 
and total number of particles (AO between 52 and 56. Note that the 
second and third spin-1 triplets are twofold degenerate. 



particle number sector N - L + \. 

Similar results can be obtained with the regular MERA. For 
instance, one can obtain an approximation to the ground state 
of a given particle number sector by adding a chemical poten- 
tial term -fi 2s to the Hamiltonian and carefully tuning 
the chemical potential term yt until the expectation value of 
the particle number N is the desired one. However, the regu- 
lar MERA cannot garantee that the states obtained in this way 
are exact eigenvectors of N. Instead, the resulting states are 
likely to have particle number fluctuations. 

Fig. 17 shows the low energy spectrum of the Heisenberg 
model H X xx for a periodic system of L — 54 sites (or 108 



spins), including the ground state and several excited states 



2. Reduction of computational costs 

The use of U(l) invariant tensors in the MERA also results 
in a reduction of computational costs. 

First, U(l) invariant tensors, when written in the canonical 
form of Eqs. 83]|86 are block diagonal and therefore require 
less storage space. Table |H| compares the number of MERA 
coefficients that need to be stored in the regular and symmet- 
ric case, for different choices of particle number assignments 
relevant to the present examples. 

Second, the computation time required to manipulate ten- 
sors is also reduced when using U(l) invariant tensors in 
the canonical form. Fig. [T~8] shows the computation time re- 
quired for one iteration of the energy minimization algorithm 
of Ref. 18 (during which all tensors in the MERA are updated 
once), as a function of the total bond dimension x- The plot 
compares the time required using regular tensors and {/(^in- 
variant tensors. For t/(l) invariant tensors, we display the 
time per iteration for three different levels of precomputation, 
as described in the appendix. The figure shows that for suffi- 
ciently large x, using U(l) invariant tensors always leads to a 
shorter time per iteration of the optimization algorithm. 

However, in the authors' reference implementation (written 
in C++ and MATLAB), using the symmetry without precom- 
putation only reduces the computational time by about a fac- 
tor of two for the largest ^ under consideration. This is due to 
the fact that maintaining the canonical form for t/(l) invariant 
tensors still imposes a significant overhead for the values of 
X considered (notice that the gap between the cost for regular 
and symmetric tensors without precomputation is still increas- 
ing as a function of^). While the magnitude of this overhead 
is necessarily dependent on factors such as programming lan- 
guage and machine architecture, more significant gains can 
be obtained by making maximum use of precomputation (giv- 
ing computation times shorter by a factor of ten or more). 
This option, however, requires a significant amount of addi- 
tional memory (see appendix), and a more convenient mid- 
dle ground can be obtained by using a partial precomputation 
scheme. 
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X = 24 



X=28 
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FIG. 18: Computation time (in seconds) for one iteration of the 
MERA energy minimization algorithm, as a function of the bond 
dimension^. For sufficiently large %, exploiting the (7(1) symmetry 
leads to reductions in computation time. The horizontal line on this 
graph shows how this reduction in computation time equates to the 
ability to evaluate MERAs with a higher bond dimension x- For the 
same cost per iteration incurred when optimising a standard MERA 
in MATLAB with bond dimension x = 20, one may choose instead 
to optimise a (7(1) symmetric MERA with partial precomputation 
and^ = 24, or with full precomputation and^ = 28. 



VI. CONCLUSIONS 

In this paper we have provided a detailed explanation of 
how a global internal Abelian symmetry may be incorporated 
into any tensor network algorithm. Following Ref. [66] we 
considered tensor networks constructed from tensors which 
are invariant under the action of the internal symmetry, and 
showed how each tensor may be decomposed according to a 
canonical form into degeneracy tensors (which contain all the 
degrees of freedom that are not affected by the symmetry) and 
structural tensors (which are completely determined by the 
symmetry). We then introduced a set of primitive operations 
P which may be used to carry out tensor network algorithms 
such as MPS, PEPS, and MERA, and showed how each of 
these operations can be implemented in a way such that the 
canonical form is both preserved and exploited for computa- 
tional gain. 

We then demonstrated the implementation of this decom- 
position for tensors with an internal t/(l) symmetry, and com- 
puted multiple benchmarks demonstrating the computational 
costs and speed-ups inherent in this approach. We found that 
although maintaining the canonical form imposed additional 
costs when combining or splitting tensor indices, for simula- 
tions of a sufficiently large scale these costs can be offset by 
the gains made when performing permutations, matrix multi- 
plications, and matrix decompositions. 

Finally, we implemented the MERA on a quantum spin 
chain with U(l) symmetry and showed that exploitation of 



this symmetry can lead to a decrease in computational cost by 
between ten and twenty times. 

To demonstrate the practical nature of these gains, we ap- 
plied U(l) symmetry to an implementation of the Multi-scale 
Entanglement Renormalization Ansatz on a quantum spin 
chain, and achieved performance increases by a factor of ten 
or more. These gains may be used either to reduce over- 
all computation time or to permit substantial increases in the 
MERA bond dimension x- 

Although in this paper we have focused on an example 
which is a continuous Abelian group, the formalism presented 
may equally well be applied to a finite Abelian group. In par- 
ticular, let us consider a cyclic group Z q , q e Z + !^J As in 
the case of t/(l), the Hilbert space decomposes under the ac- 
tion of the group into a direct sum of one dimensional irreps 
which are each characterized by an integer charge a, and con- 
sequently most of the analysis presented in this paper remains 
unchanged. In particular, matrices which are invariant under 
the action of the group will be block diagonal in the basis la- 



beled by the charge a according to Eq. 35 and symmetric ten- 
sors enjoy the canonical decomposition stated in Eqs. 85p6 
The only objects which need modification are the fusion and 
splitting maps, which need to be altered so that they encode 
the fusion rules of Z q instead of U(\). For a cyclic group Z q , 
the fusion of two charges a and a' gives rise to a charge a" ac- 
cording to a" = {a + a')\ q where L indicates that the addition is 
performed modulo q. For example, Z3 has charges a — 0, 1,2, 
and the fusion rules for Z3 take the form a x a' — > a" where 
the value of a" is given in the following table: 







a 









1 


2 








1 


2 


a' 1 


1 


2 





2 


2 





1 



More generally, a generic abelian group will be charac- 
terised by a set of charges (01,02,03, . . .). When fusing two 
such sets of charges (01,02,03, . . .) and {a\,a' 2 ,a' v . . .), each 
charge a,- is combined with its counterpart a' i according to the 
fusion rule of the relevant subgroup. Once again, this be- 
haviour may be encoded in a single fusion map T lllsc and its 
inverse Y spl \ The formalism presented in this paper is there- 
fore directly applicable to any Abelian group. 
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APPENDIX: USE OF PRECOMPUTATION IN ITERATIVE 
ALGORITHMS 

We have seen that the use of the canonical form given in 
Eqs. 85p6 to represent t/(l) invariant tensors can potentially 
lead to substantial reductions in memory requirements and in 
calculation time. We also pointed out, however, that there is an 
additional cost in maintaining an invariant tensor in its canoni- 
cal form, and that this is associated with the reshaping (fusing 
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and/or splitting) of its indices. In some situations this addi- 
tional cost may significantly reduce, or even offset, the bene- 
fits of using the canonical form. 

In this appendix we investigate techniques for reducing this 
additional cost in the context of iterative tensor network algo- 
rithms. Many of the algorithms discussed in Sec. HE are it- 
erative algorithms, repeating the same sequence of tensor net- 
work manipulations many times over. Examples include algo- 
rithms which compute tensor network approximations to the 
ground state by minimizing the expectation value of the en- 
ergy, or by simulating evolution in imaginary time, with each 
iteration yielding an increasingly accurate approximation to 
the ground state of the system. 

The goal of this appendix is to identify calculations which 
depend only on the symmetry group, and are independent of 
the variational coefficients of such algorithms. Where these 
calculations are repeated in each iteration of the algorithm, 
we can effectively eliminate the associated computational cost 
by performing them only once, either during or prior to the 
first iteration of the algorithm, and then storing and reusing 
these precomputed results in subsequent iterations. We will 
illustrate this procedure by considering the precomputation of 
a series of operations applied to a single tensor T. 

To do this, we begin by revisiting the fusion and splitting ta- 
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bles of Sec. Ill C and introducing a graphical representation of 
these objects. We then introduce a convenient decomposition 
of a symmetric tensor into a matrix accompanied by multiple 
fusion and/or splitting tensors, and linear maps F that map one 
such decomposition into another. These linear maps are inde- 
pendent of the coefficients of the tensor being reorganized, 
and consequently they are precisely the objects which can be 
precomputed in order to quicken an iterative algorithm at the 
expense of additional memory cost. Finally we describe two 
specific precomputation schemes, differing in what is precom- 
puted and in how the precomputed data is utilized during the 
execution of the algorithm, in order to illustrate the trade off 
between the amount of memory needed to store the precom- 
putation data and the associated computational speedup which 
may be obtained. In practice, the nature of the specific imple- 
mentation employed will depend on available computational 
resources. 



1. Diagrammatic notation of fusing and splitting tensors 

In describing how we can precompute repeated manipula- 
tions of this tensor T, we will find it useful to employ diagram- 
matic representations of the fusion and splitting tables T fuse 



and Y spli ' introduced in Sec. IIIC These tables implement a 
linear map between a pair of indices and their fusion product, 
and thus can be understood as trivalent tensors having two in- 
put legs and one output leg (or vice versa) in accordance with 
Sec. |II F| We choose to represent them graphically as shown in 
Fig.|19|i), where the arrow within the circle always points to- 
ward the coupled index. The linear maps T rilSL and Y spli ' are uni- 



tary, and consequently we impose that the tensors of Fig. 19 i) 



must satisfy the identities given in Fig. 19 'ii), corresponding 
to unitarity under the action of the conjugation operation erri- 
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FIG. 19: (i) Graphical representation of the fusion tensor T f " SE and 
the splitting tensor T spl ". (ii) The tensors T fuse and T spl " are unitary, 
and thus yield the identity when contracted pairwise as shown, (iii) 
A fusion tensor decomposed into two parts. The first part (indicated 
by a circle with an arrow) performs the tensor product of input ir- 
reps, n A tA x n B t B . The result is an index that labels pairs (iiAt A , n B t B ). 
The second part (indicated by a rectangle) is a permutation that as- 
sociates each pair (n A tA, n B t B ) with a unique (tiA B t„ AB ), corresponding 
to a vector of the coupled basis of 15 



ployed in diagrammatic tensor network notation (vertical re- 
flection of a tensor and the complex conjugation of its com- 
ponents, typically denoted '). Our notation also reflects the 
property, first noted in section IIIC that T ruse and T split may 



be decomposed into two pieces (Fig. 19 lii)). For the fusion 
tensor, we identify the first piece (represented by a circle con- 
taining an arrow) with the creation of a composed index using 
the manner we would employ in the absence of symmetry Q. 
The second piece, represented by the small square, permutes 
the basis elements of the composed index, reorganizing them 
according to total particle number. The two components of the 
splitting tensor are then uniquely defined by consistency with 
the process of conjugation for the diagrammatic representa- 



tion of tensors, and with the unitarity condition of Fig. 19 ii) 



These requirements have an important consequence. Sup- 
pose the first part of T fl,sc implements b x c — > d by iterating 
rapidly over the values of b and more slowly over the values 
of c, and b lies clockwise of c on the graphical representation 
of T fusc . This then means that on the graphical representa- 
tion of T spli ' which implements d — > b x c, index b must lie 
coM«ferclockwise of c. It is therefore vitally important to dis- 
tinguish between the splitting tensor and a rotated depiction 
of the fusing tensor. To this end we require that when us- 
ing this diagrammatic notation, all tensors (with the exception 
of the fusion and splitting tensors) must be drawn with only 



downward-going legs, as seen for example in Fig. 20 though 
the legs are still free to carry either incoming or outgoing ar- 
rows as before. 
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FIG. 20: Binary tree decomposition of a symmetric tensor f having 
components T jlhj}i4 j 5 j 6 . The tree T is comprised of a matrix M as 
the root node, four splitting tensors as internal nodes, and if, 
as its leaf indices. No incoming or outgoing arrows are indicated on 
the indices in the figure, as the decomposition is valid for any such 
assignment of directional arrows. 





h h h h 



FIG. 21: Two possible tree decompositions of a rank-4 tensor f. 
Different choices T\,Tz,--- of tree decomposition for tensor T lead 
to different matrices M { ,Mi, ■ ■ ■ for the same tensor. 
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FIG. 22: Tree decompositions of tensor f are obtained by contracting 
the tensor with an appropriate resolution of the identity on its indices, 
selected according to the desired choice of the fusion tree T. In each 
instance, evaluation of the contents of the shaded region yields the 
appropriate matrix M. 



2. Tree decomposition 

We find it convenient to decompose a rank-£, i/(l) invari- 
ant tensor T, having components 7 1 ,, as a binary tree ten- 
sor network T consisting of a matrix M which we will call 
the root node, and of k - 2 splitting tensors T splil as branching 
internal nodes, with the leaf indices of tree T corresponding 
to the indices \i\,h, ■ ■ ■ , ik) of tensor T. We refer to decom- 
position T as a tree decomposition of T . 



Fig. 
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shows an 



example of tree decomposition for a rank-6 tensor. It is of the 
form 



T 

H'2 



- z 



M hi :\' : ' . . Y spli ' . . T spl " 

J' J 2 h->H,J3 j2->74,i6 Jl-> 



(102) 

where {j\,j2, 73, ]a) are the internal indices of the tree. 

The same tensor T may be decomposed as a tree in many 
different ways, corresponding to different choices of the fu- 
sion tree. As an example we show some different, but equiv- 
alent, decompositions of a rank-4 tensor in Fig. 21 Differ- 
ent choices T\, 7~2, • • • of tree decomposition for tensor T will 
lead to different matrices representations M\,M^, ■ ■ • of the 



same tensor. Finally, Fig. 22 shows how to obtain the tree de- 
compositions from r (| ,' 2 ,',,' 4 by introducing an appropriate res- 
olution of the identity, constructed from pairs of fusion op- 
erators Y 1 "™ and splitting operators T splil in accordance with 



The representation of a tensor f by means of a tree decom- 
position is particularly useful because many tensor network 
algorithms may be understood as a sequence of operations 
carried out on tensors reduced to matrix form. For example, 
tensor network algorithms such as MPS, MERA, and PEPS 
consist primarily of (i) tensor network contractions, and (ii) 
tensor decompositions. In Sec. II D we argued that all such 



operations may be reduced to matrix multiplications, matrix 
decompositions, and a set of primitive operations P. When 
tensors are updated in these algorithms they are typically cre- 
ated as matrices, to which operations from P are then applied, 
and when they are decomposed or contracted with other ten- 
sors, this once again may take place with the tensor in matrix 
form. Any such matrix form may always be understood as 
the matrix component of an appropriate tree decomposition T 
of tensor T, where the sequence of operations required to re- 
shape tensor T to matrix M corresponds to the contents of the 



shaded area in Fig. 22 



Fig. 19 ;ii). 



3. Mapping between tree decompositions 

Suppose now that we have a tensor T in matrix form M\, 
which is associated with a particular choice of tree decompo- 
sition 7~i , and we wish to transform it into another matrix form 
M2, corresponding to another tree decomposition To- As indi- 
cated, this process may frequently arise during the application 
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M, 



FIG. 23: A matrix Mj can be reorganized into another matrix M 2 
by means of fusion tensors, splitting tensors, and the permutation of 
indices. These operations define a one to one linear map T that acts 
to reorganize the coefficients of Mi . T does not depend on the coef- 
ficients of Mi, but solely on the sequence of operations performed. 



of many common tensor network algorithms. The new matrix 
M2 can be obtained from M\ by means of a series of reshap- 
ing (splitting/fusing) and permuting operations, as indicated 
in Fig. [23] and this series of operations may be understood as 
defining a map T: 



M 2 = r(M0. 



(103) 



The map F is a linear map which depends only on the tree 
structure of T\ and T2, and is independent of the coefficients 
of M\ . Moreover T is unitary, and it follows from the con- 
struction of fusing and splitting tensors and the behaviour of 
permutation of indices (which serves to relocate the coeffi- 
cients of a tensor) that F simply reorganizes the coefficients of 
M\ into the coefficients of M 2 m a one-to-one fashion. 

Therefore, one way to compute the matrix M 2 from matrix 
Mi is by first computing the linear map F, which is indepen- 
dent of the specific coefficients in tensor T, and by then ap- 
plying it to Mi . 



4. Precomputation schemes for iterative tensor network 
algorithms 

The observation that the map F is independent of the spe- 
cific coefficients in M; is particularly useful in the context of 
iterative tensor network algorithms. It implies that, although 
the coefficients in Mi will change from iteration to iteration, 
the linear map F in Eq. |103| remains unchanged. It is therefore 
possible to calculate the map F once, during the first iteration 
of the simulation, and then to store it in memory and re-use it 
during subsequent iterations. We refer to such a strategy as a 



precomputation scheme. Fig. 24 contrasts the program flow of 
a generic iterative tensor network algorithm with and without 
precomputation of the transformations F. 

Using such a precomputation scheme a significant speed-up 
of simulations can be obtained, at the price of storing poten- 
tially large amounts of precomputed data (as a single iteration 
of the algorithm may require the application of many different 
transformations F). There therefore necessarily exists a trade- 
off between the amount of speed-up which can be obtained 
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k = count^— »— ( End ) 
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FIG. 24: Flow diagram for the execution of a predetermined num- 
ber of iterations of a generic iterative tensor network algorithm (i) 
without any precomputation and (ii) with precomputation of the op- 
erations r. 



and the memory requirement which this entails. In this sec- 
tion we describe two different precomputation schemes. The 
first one fully precomputes and stores all maps F, and is rela- 
tively straightforward to implement. This results in the maxi- 
mal increase in simulation speed, but implementation requires 
a large amount of memory. The second scheme only partially 
precomputes the maps F, resulting in a moderate speed-up of 
simulations, but with memory requirements which are also 
similarly more modest. 



a. Maximal precomputation scheme 

As noted in Sec.|3]of this appendix, applying the map F to 
a matrix M\ simply reorganizes its coefficients to produce the 
matrix M 2 . Moreover, if the indices of matrices M\ and M 2 
are fused to yield vectors V\ and V2 then the map F may be 
understood as a permutation matrix, and this in turn may be 
concisely represented as a string of integers F = y\, . . . , J\m x \ 
such that entry i of V 2 = FV\ is given by entry y -, of vector 
V*! . Because all of the elements from which F is composed are 
sparse, unitary, and composed entirely of zeros and ones, the 
permutation to which F corresponds may be calculated at a to- 
tal cost of only <9(|Mi|), where \M\ \ counts only the elements 
of Mi which are not fixed to be zero by the symmetry con- 
In essence, 



straints of Eq. 85 In essence, for each element of the vector 
Vi one identifies the corresponding number and degeneracy 
indices {nf x , tf' ) on each leg 2 € {1, 2} of matrix M\ . One can 
then read down the figure, applying each table T fl,sc or T spl " 
in turn to identify the corresponding labels («',f') on the in- 
termediate legs, until finally the corresponding labels on the 
indices of M 2 are obtained. There is then a further 1:1 map- 
ping from each set of labels (nf 2 ,tf 2 ), (nf 2 , t™ 2 ) on M 2 to the 
corresponding entry in V 2 , completing the definition of F as a 
map from V\ to V 2 . 

Storing the map F for a transformation such as the one 
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FIG. 25: (i) Permutations applied to one or more legs of a fusion 
or splitting tensor can be replaced by an appropriate permutation on 
the coupled index. This process can be used to replace all permuta- 
tions applied on internal indices of a diagram such as Fig. |23|with 
net permutations on the indices of M\ and on the open indices of 
the network, as in shown in (ii). The residual fusion and splitting 
operations, depicted as an arrow in a circle, simply perform the ba- 
sic tensor product operation and its inverse, {2}-{3]l as described in 
Fig.fMiii). 



shown in Fig. 23 imposes a memory cost of 0(\M\\). The 
application of this map also incurs a computational cost of 
0(\M\\), but computational overhead is saved in not having to 
reconstruct the map F on every iteration of the algorithm. 



b. Partial precomputation scheme 

The 0(\M\\) memory cost incurred in the previous scheme 
can be significant for large matrices. However, we may re- 
duce this cost by replacing the single permutation Y employed 
in that scheme with multiple smaller operations which may 
also be precomputed. In this approach M\ is retained in matrix 
form rather than being reshaped into a vector, and we precom- 
pute permutations to be performed on its rows and columns. 

First, we decompose all the the fusion and splitting ten- 
sors into two pieces in accordance with Fig. 



19 iii). Next, 



a single permutation applied to the coupled index (Fig. 25 i)). 
We use this to replace all permutations on the intermediate in- 
dices of the diagram with equivalent permutations acting only 
on the indices of M\ and the open indices, as shown for a sim- 
ple example in Fig. 25 ii). The residual fusion and splitting 



operations, depicted by just a circle enclosing an arrow, then 
simply carry out fusion and splitting of indices as would be 
performed in the absence of symmetry Q-([3j. These opera- 
tions are typically far faster than their symmetric counterparts 
as they do not need to sort the entries of their output indices 
according to particle number. 

In subsequent iterations, the matrix Mi is obtained from M\ 
by consecutively 

1 . Permuting the rows and columns of M\ using the pre- 
computed net permutations which act on the legs of M\ . 

2. Performing any elementary (non-symmetric) splitting, 
permuting of indices, and fusing operations, as de- 
scribed by the grey-shaded region in Fig.|25jn). 

3. Permuting the rows and columns of the resulting matrix, 
using the precomputed net permutations which act on 



the open legs of Fig. 25 li) 



When matrix Mi is defined compactly, as in ( 85 i, so that el- 
ements which are identically zero by symmetry are not ex- 
plicitly stored, a tensor T is constructed from multiple blocks 
identified by t/(l) charge labels on their indices (f n)m ... nk in 



we recognise that any permutations applied to one or more 
legs of a fusion or splitting tensor may always be written as 



Eq. 85 1. Under these conditions the elementary splitting, fus- 
ing and permutation operations of step 2 above are applied 
to each individual block, but some additional computational 
overhead is incurred in determining the necessary rearrange- 
ments of these blocks arising out of the actions performed. 
This rearrangement may be computed on the fly, or may also 
be precomputed as a mapping between the arrangement of 
blocks in M\ and that in M 2 . 

The memory required to store the precomputation data in 
this scheme is dominated by the size of the net permutations 
collected on the matrix indices, and is therefore of 0( V|Mi|). 
The overall cost of obtaining Mi from Mi is once again of 
0(\M\\), but is in general higher than the previous scheme 
as this cost now involves two complete permutations of the 
matrix coefficients, as well as a reorganisation of the block 
structure of M\ which may possibly be computed at runtime. 
Nevertheless, in situations where memory constraints are sig- 
nificant, partial precomputation schemes of this sort may be 
preferred. 
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